Gravitational wave asteroseismology with fast rotating neutron stars 
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We investigate damping and growth times of the quadrupolar /-mode for rapidly rotating stars 
and a variety of different polytropic equations of state in the Cowling approximation. This is the 
first study of the damping /growth time of this type of oscillations for fast rotating neutron stars in a 
relativistic treatment where the spacetime degrees of freedom of the perturbations are neglected. We 
use these frequencies and damping/growth times to create robust empirical formulae which can be 
used for gravitational wave asteroseismology. The estimation of the damping/growth time is based 
on the quadrupole formula and our results agree very well with Newtonian ones in the appropriate 
limit. 
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I. INTRODUCTION 

During the birth of a proto-neutron star or the merg- 
ing of two older compact stars, violent non-radial oscil- 
lations may be excited, resulting in the emission of sig- 
nificant amounts of gravitational radiation [1]. The de- 
tection of gravitational waves from oscillating neutron 
stars will allow the study of their interior, in the same 
way as helioseismology provides information about the 
interior of the Sun. It is expected that the identification 
of specific pulsation frequencies in the observational 
data will reveal the true properties of matter at densities 
that cannot be probed today by any other experiment. 
In this paper, we present new empirical relationships for 
mode-frequencies and damping times of the quadrupo- 
lar /-mode for rapidly rotating neutron stars, extending 
previous studies which deal with the non-rotating case 

nail. 

These original suggestions about gravitational wave 
asteroseismology have been supported by many com- 
plementary works which studied specific features of os- 
cillation spectra for various compact objects, such as 
t5^ical neutron stars |,5-;12|, but also for strange IITSHTSi 
or superfluid stars 1161 [17|. More recently, it has also 
been suggested that one may use asteroseismology to 
find the imprints of scalar or even vector components 
of gravity IITSHZTJI . It should be noted here that all previ- 
ous studies IITI- ET] have been performed for non-rotating 
relativistic stars. The treatment of rotation was always 
a problem in general relativity and thus the majority 
of the studies for the oscillation spectra of fast-rotating 
compact stars was done mainly in Newtonian theory 
w^hich gives only qualitative answers. 

Since stellar oscillations may become unstable in the 
presence of rotation, there was an increased interest dur- 
ing the last decade or so to study the dynamics of rotat- 
ing stars, also thanks to the discovery of the r- and w- 
mode instability Ii22i -t26 1 . Still, the majority of these stud- 



ies have been performed in Newtonian theory Il27f|29l 
while there are only a few works in which GR has been 
used; mainly in the so called slow-rotation approximation. 
The slow-rotation approximation was successfully ap- 
plied to study various aspects of the r-mode instabil- 
ity Il30l |3TI , effects of uniform and differential rotation 
on the oscillation spectrum [i32n34J and on the crustal 
modes 1 



As it has been originally suggested by Chandrasekhar 
1 36 1 and verified by Friedman & Schutz |37, 38| certain 
non-axisymmetric pulsation modes may grow exponen- 
tially in rotating stars; this is due to the emission of grav- 
itational waves and is called CFS instability. Exploring 
this type of instability in rapidly rotating stars turned 
out to be very difficult. In linear perturbation theory 
for example, rapid rotation was never treated properly 
until recently; almost all formulations of the relevant 
perturbation equations were prone to numerical insta- 
bilities either at the surface or along the rotation axis 
of the neutron star. Thus, it was not surpising that the 
first results for the oscillations of rapidly rotating stars 
were derived using evolutions of the non-linear equa- 
tions [39 -42 J. Still all these studies were purely axisym- 
metric and thus the effects of rotation on the spectra was 
present only for very high rotation rates. Rotational in- 
stabilities are driven by non-axisymmetric modes and 
thus these first 2D calculations where not of much use 
for their study. 

In the last two years there was significant progress in 
the study of non-axisymmetric perturbations of rapidly 
rotating neutron stars. For the first time it was possible 
to calculate in GR the oscillation spectra of fast rotating 
relativistic stars by using the linearized form of the fluid 
equations. Thus the effect of fast rotation on /- and r- 
modes has been demonstrated while the critical points 
for the onset of the /-mode (CFS) instability have been 
derived 1431 . In addition it has been demonstrated that 
there is a way to derive empirical relations connecting 



the oscillation frequencies with the rotation of the stars. 
This study has been recently extended to ^-modes |44J 
and even more recently has been expanded to study the 
oscillation spectra of fast and differentially rotating neu- 
tron stars |45|. 

It should be noted that the previous results have 
been derived using the so-called Cowling approxima- 
tion where the spacetime is assumed to be frozen. This 
approximation is very good for the estimation of the 
spectra of r- and g-modes but it gives only qualitatively 
good results for the /-mode. Moreover, using non-linear 
codes it became possible for the first time to study the 
complete problem |46|, i. e. the non-axisymmetric stellar 
oscillations of fast rotating stars without the constraints 
of the Cowling approximation. The results are in qual- 
itative agreement with those found in Il43l and for the 
critical point for the onset of the /-mode instability with 
the studies presented in [47 1. 

The next step for gravitational wave asteroseismol- 
ogy is to use additional information about the damping 
times to construct model-independent relations which 
allow for a robust determination of stellar key param- 
eters. The damping time t of the potentially CFS- 
unstable branch in the high rotation regime for example 
can be approximated very accurately by using a sim- 
ple relation of the form tq/t = sgn((j,)0.256(cr,/(Jo)^, 
w^here tq and (jg are the damping time and mode fre- 
quency of the nonrotating model respectively, sgn (x) is 
the signum function and u, is the actual mode frequency 
in the inertial frame. 

The structure of the paper is as follows. In Section 
In] we give an essential overview about our method of 
computing mode-frequencies and damping times of the 
/-mode. We then show the results of our simulations in 



Section III where we present empirical relations which 
can potentially be used to estimate masses, radii and 
angular frequencies of rapidly rotating neutron stars. 
A more elaborate discussion about the numerical pro- 
cedure, the equations of state and background models 
used in this study as well as a consistency check can be 
found in the Appendix. 



II. PROBLEM SETUP 

Mode frequencies and damping times of neutron star 
oscillations can be calculated in two different ways. In 
a time-independent boundary-value formulation of the 
perturbation equations, they are directly obtained from 
the real and the imaginary part of the complex eigen- 
frequencies. In a time-dependent evolution problem on 
the other hand, both these quantities have to be com- 
puted in post-processing routines. The frequency of spe- 
cific fluid modes is obtained by Fourier-transforming 
the time-series at different points inside the star into the 



frequency domain and correctly identifying the corre- 
sponding peaks in the power spectrum, see |43| for a 
detailed description of this method. 

Concerning the damping time, one has to calculate 
both the energy that is contained within a mode as 
well as the energy loss due to gravitational wave emis- 
sion w^hich in our case is done via the quadrupole for- 
mula, see e. g. |28| for an application of this formal- 
ism to r-modes. A different procedure is to use the be- 
haviour of metric perturbations at future null infinity to 
derive a gauge-invariant expression for the luminosity, 
see |48J. However, since we will work in the Cowling- 
approximation, the quadrupole formalism is utilized in 
this study. 

The energy of a mode in a comoving frame is then 
given by 
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where p is the rest-mass density and 6p its correspond- 
ing perturbation, 6p, 5^ the perturbations of pressure 
and gravitational potential respectively. On the other 
hand, the quadrupole formula for the emission of grav- 
itational radiation gives 
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is the coupling constant for spherical mode number /, ct, 
the mode frequency in the inertial frame, n the angular 
velocity of the neutron star and where SDi,„, S}]„, are the 
mass- and the current-multipole moments respectively. 
The damping time is then computed by 
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Several remarks are now worth considering. First, as 
already mentioned we are working in the Cowling- 
approximation w^hich means that the spacetime is kept 
fixed during the time-evolution, i. e. <50 = in equa- 
tion |ll|. Second, we will focus on the nonaxisymmet- 
ric / = \m\ = 2 fundamental mode since in general 
it has the smallest growth-time and is therefore more 
relevant in real astrophysical scenarios. Higher order 
modes typically not only grow on larger timescales but 
are also damped stronger by various dissipative effects. 
In the case of pressure modes, the emission of grav- 
itational waves is to a great extent due to the mass- 
quadrupole moment and we will neglect the current 
quadrupole moment S}22- Taking these comments into 



account, equations ||T| and ||2} adjusted for /-modes in 
the Cowling-approximation, read 
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as mass-quadrupole moment. 

The damping time t^,- depends crucially on how the 
mode frequencies of non-axisymmetric perturbations 
change with rotation rate. While degenerated in the 
non-rotating limit, the frequencies of modes with the 
same spherical mode number / but opposite azimuthal 
index m = ±\m\,i. e. co- and counterrotatrng modes, di- 
verge. Figure [1] shows an example of this behaviour for 
a certain sequence of equilibrium models with increas- 
ing angular frequency n/27r in a system comoving with 
the star. The power spectral density of the pressure per- 
turbation variable H is taken at an arbitrary point inside 
the neutron star (typically s = t = 0.5; for a description 
of the computational domain and the coordinates used 
there, see Appendix IA| and colour-coded for the differ- 
ent models. In the nonrotating limit, one can identify 
various peaks with the strongest and sharpest ones lo- 
cated at cri/2Ti = 3.837 kHz and at i72/2n = 9.432 kHz. 
Inspection of the corresponding eigenfunctions shows 
that the peak at Vi belongs to the quadrupolar /-mode 
while 1/2 matches its first overtone, the ^pi-mode with 
an additional node in radial direction. Increasing the 
angular velocity leads to a splitting of nonaxisymmetric 
perturbations which can most clearly be seen for the two 
modes mentioned. Other peaks in Figure IT] split as well, 
however due to broader edges this behaviour is harder 
to follow. 

The imaginary part of the complex mode frequency, 
which is approximated by equation ||4|, controls the 
exponential damping or growing of non-axisymmetric 
perturbations. For nonrotating stars it is E > and 
dE/dt < 0, see equations Q and l|2}. The imaginary 
part of the mode frequency then becomes positive, indi- 
cating a damped oscillation. The perturbations remain 
damped as long as ct-j/ct, + mCi) > 0, that is counter- 
rotating modes in the comoving frame are still counter- 
rotating in the inertial frame. This behaviour changes 
once the pattern speed of the mode is matched by the 
angular velocity. In this case, dE/dt = and there is no 
loss of energy due to gravitational radiation. Finally, for 
background configurations that allow mode frequencies 
beyond the zero-frequency limit in the inertial frame it is 




'^ (kHg) 



FIG. 1. Splitting of the power spectral density (normalized 
units) for non-axisymnietric \m\ = 2-niodes in a comoving 
reference frame with corresponding mode-frequency u^. The 
mass-shedding limit for this particular sequence is reached for 
n/2n = 2.18 kHz. 



dE/dt > and the oscillation is exponentially growing 
on a timescale given by Tgr- 

These are the astrophysically most interesting cases 
since the oscillation is unstable in this regime, emitting 
significant amounts of gravitational radiation. While 
damping times have already been computed for non- 
rotating stars fSj-^ IZl Il9 1 as well as in the relativistic 
slow-rotation approximation |50|, so far there are no 
numerical simulations for calculating damping times of 
rapidly rotating relativistic models. However, in New- 
tonian theory the quadrupole formula has been used 
successfully for computing damping times of rotating 
polytropes IISTl |52| and recently also for gravitational- 
wave extraction of rotating superfluid stars 153J. 



III. ASTEROSEISMOLOGY 



A. Frequencies of Co- and Counterrotating Modes 



In 1431 we already presented results about the split- 
ting of the fundamental mode in a coordinate frame co- 
moving with the star. The conclusions there suggested, 
that while in the inertial frame the exact details of the 
/-mode splitting depend on the particular background 
model and the given equation of state, when properly 
normalized in the comoving frame, one can actually 
construct a model-independent relation between mode 
frequency and rotation rate. 
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FIG. 2. Co- and counterrotating branches of the different polytropic EoS studied in this paper. Panel (i): Mode frequencies in 
the inertial frame; the solid line of each EoS depicts the more compact sequence while the dashed line traces the less compact 
configurations. Panel (ii): Normalized mode frequencies and fitting curve in the comoving frame; the larger circles represent 
again the more compact models while the small circles stand for the less compact ones, og is the frequency in the nonrotating 
Umit, n^ represents the Kepler-limit. 



In parts, this is due to the following reason. In the in- 
ertial frame there is a clear cut between configurations 
that allow for potentially CFS-unstable models and con- 
figurations which never become CFS-unstable. If one 
were to fit for example all counter-rotating modes with 
just one fitting curve, either all models will become un- 
stable at some point or no model at all. Clearly, this 
would be implausible for a proper fitting. 

In the comoving frame on the other hand, the situa- 
tion is different. There, the two branches of counter- and 
corotating modes never reach the zero-frequency limit 
even remotely and this is due to the fact that the abso- 
lute value of the splitting between the m = \m\ and the 
m = —\m\ branches is smaller there. The relation be- 
tween mode frequencies u, in the inertial frame and a^ 
in the comoving frame is given by 



otO. 



(8) 



It then follows from equation ||8), that the frequency sep- 
aration A between the two branches transforms accord- 
ing to 



A; = Ac + 2\m\Ci, 



(9) 



that is the separation is always smaller in the comoving 
system and it is actually the corotating branch that has 
the lower frequencies in this frame, see e. g. [43, 54]. 

Figure [2] shows the results of our simulations regard- 
ing the mode frequencies and as expected, they show 
a large variety in the inertial frame. Depending on the 



actual configuration and equation of state, some mod- 
els become unstable before reaching the mass-shedding 
limit, some are only marginally unstable and some re- 
main stable even at the Kepler frequency. 

Despite this apparent diversity in the inertial frame, 
the two branches can very well be fitted with a second 
order polynomial in the comoving frame, see Panel (//) 
of Figure [2] There, the mode frequency is normalized 
by its value in the nonrotating limit tTg while the angular 
velocity is prescribed in units of the Kepler frequency 
rjj^. As one can see, the fitting describes the overall be- 
haviour of the mode frequencies very well; only for val- 
ues close to the mass-shedding limit the various curves 
for the different EoS show larger deviations from the 
quadratic fit. One should also keep in mind that the 
order of co- and counterrotating modes is reversed in 
the comoving frame, that is while unstable modes have 
lower frequencies in the inertial frame, they represent 
the high-frequency branch in the comoving frame and 
vice versa. 

Based on the data of our simulations, we propose the 
following relationships for the quadratic fitting polyno- 
mials. It is 
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for the always stable {m = — 2) and 
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(10) 



(11) 



for the potentially unstable {m = 2) branch. This also 



agrees very well with our previous findings for a more 
limited set of equilibrium configurations and equations 
of state in 1 43 1 . 

In addition, an auxiliary condition is needed that con- 
nects the mode frequency in the nonrotating limit ug 
with fundamental stellar parameters. It is well known, 
that for the /-mode ctq scales with the mean density 
{M/R^y^^, see e. g. ||3l [A], where fitting coefficients 
based on a variety of realistic EoS are provided. 

We repeated this calculation in the time-domain and 
with our set of equation of states and the results are de- 
picted in Figure Is] In order to better compare them with 
the findings in |3|, we picked out the very same range 
of mean densities. Depending on the particular EoS, it 
may happen that certain configurations never reach the 
range of mean densities depicted in Figure pi for exam- 
ple the very stiff EoS PI. 4 is completely absent. 
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FIG. 3. Scaling of the /-mode frequency with mean density in 
the nonrotating limit. M and K are given in km, uq in kHz. 



The frequencies computed with our approach generally 
tend towards larger values. This can very well be under- 
stood with the Cowling-approximation w^hich is used 
to simplify the time-evolution equations, see Appendix 
[A| The freezing of all metric perturbations systemati- 
cally overestimates pressure-mode frequencies though 
this effect becomes less pronounced for higher order 
modes, see also [46 1 and references therein. Another 
difference is the use of polytropic equations of state in 
contrast to tabulated EoS utilized in |3|; this also has an 
effect on the mode frequencies and preliminary studies 
show that this might be an even stronger restriction than 
the Cowling-approximation ||55] . 

Nevertheless, it is still possible to fit the frequencies 
very well with a linear dependence in the mean density 



and we find 
1 



„ (Tq (kHz) = 0.498 -F 2.418 , -, 



M 



1/2 



where we introduced the dimensionless variables 
M , ^ R 



M 



1.4 Mrr 



and R 







10 km 



(12) 



(13) 



B. Damping Times of Co- and Counterrotating Modes 

A similar procedure can be applied to the damping 
times of the two fundamental mode branches. However, 
in order to find a model-independent relation, we can- 
not use the angular velocity O directly as a measure of 
the rotation rate as it was done in Panel (ii) of Figure 
The reason for this is, that the damping time changes 
its sign w^hen a particular fundamental mode eventu- 
ally becomes unstable, see the discussion in Section In] 
A negative damping time signals an exponential growth 
instead of a damped oscillation. If one were to fit the 
various damping times of the counterrotating modes as 
function of angular velocity with just one fitting curve, 
then again either all models would become unstable at 
some point or no model at all; see also the discussion 
in Section IIII Al w^here we discussed a similar effect for 
mode frequencies in the rnertial frame. 

What is needed for the potentially unstable branch 
is a quantity that also changes its sign when a mode 
becomes prone to the CFS-instability and which is a 
monotonic function of the rotation rate. The /-mode fre- 
quency in the inertial frame cr, exactly conforms to these 
requirements. From equations iJS), l]6) we can also make 
an estimation on how the damping times depend on the 
mode frequency cr. Since dE/dt ^ a^ and for any oscil- 
lation £ ^ (7^, we have 



1 dE/dt 



rsj (7 



(14) 



In Panel (f) of Figure H] we show the corresponding re- 
sults of our simulations. There we plot normalized val- 
ues of (1/t)^'* against normalized mode frequencies (7/ 
in the inertial frame. Here, the normalization constants 
are the damping time tq and mode frequency (Jq in the 
nonrotating limit. The rotation rate increases from right 
to left, i. e. with decreasing frequency, where the points 
Vi = (1, 1) and 1^2 = (0,0) correspond to the nonrotat- 
ing case and to an infinite damping time at the onset of 
the CFS-instability respectively. 

Since for a linear fitting, the two fixed points Vi and 
p2 already determine the fitting coefficients indepen- 
dent of our actual simulations, we decided to fit the 
data points with the next highest reasonbable polyno- 
mial w^hich would be of cubic order This not only fits 
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FIG. 4. Model-independent relations for the damping time. Larger circles represent the more compact models while small circles 
stand for the less compact configurations of each EoS. Panel (i): Damping times for the counterrotating branch, a" is the mode 
frequency in the inertial frame. Panel (ii): Damping times for the corotating branch, a^ is the mode frequency in the comoving 
frame. 



the simulation data better but using a third order poly- 
nomial also turns out to be a very good approximation 
for the damping times of the corotating branches as we 
will see later 

A generic cubic pol5n-iomial is of the form i/{x) = 
ax^ + bx^ + ex + d. Imposing the constraints, that the 
fit has to pass through Vi and V2 leads to d = and 
c = 1 — a — h and least-square methods finally yield 
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w^here sgn (x) is the sign function. Equations lIS} and 
| (Tl) can then be used to cast this relation into a form 
that depends on the angular velocity again. 

For the corotating branch, this type of scaling will not 
work for the following reasons. First, the damping times 
of the stable branch decrease with the rotation rate so 
instead of fitting 1/t we will rather use t itself. Sec- 
ond and more importantly, the frequencies of the coro- 
tating branch do not depend monotonically on the rota- 
tion rate, see Panel (//) of Figured] For angular veloc- 
ities near the mass-shedding limit, the mode frequen- 
cies tend to decrease. However, as already discussed 
in Section III A the order of the two nonaxisymmetric 
branches is reversed in the comoving frame. There, the 
frequencies of the comoving modes indeed again de- 
crease monotonically with rotation rate, see once more 



Panel {ii) of Figure|2] Thus, for these modes we can use 
the mode frequencies in the comoving frame as indica- 
tor for the angular velocity and the results for this type 
of parametrization is depicted in Panel [ii) of Figure El 
Again, the rotation rate increases from right to left with 
Vi = (1,1) representing the nonrotating limit. 

Here, the spread of the data points in the [a, T)-plane 
is larger when compared to the unstable branch but still 
it can be fitted very well with a third order pol3n-iomial. 
Especially the boosted decrease in the damping times 
for high rotation rates w^hich directly correlates with the 
decrease of the mode frequencies in the inertial frame is 
captured very good with a cubic fit and cannot be repro- 
duced properly by a quadratic poynomial. 

Starting again with a generic cubic fitting function 
and including the point Vi leads to 
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(16) 



where equations l|8) and ( [10^ can be used to replace the 
comoving mode frequency by the rotation rate of the 
star. 

Similar to Section |III Aj an additional, model- 
independent relation for the damping time Tq of non- 
rotating configurations is needed. As it was shown in 
[3 1, the behaviour of K^/ (M'^t) with respect to the com- 



pactness M/R proves to be quite insensitive to details 
of the particular equation of state. 
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FIG. 5. Scaling of the /-mode damping time in the nonrotating 
Limit. M, R and tq are given in km. 

In Figure [5] we again compare our results with the cor- 
responding findings in [3J. In this case, only the soft- 
est equations of state from our sample, i. e. EoS II, A, 
P0.66 and Pl.O attained high values in M/R which al- 
low for enough data points and a proper fitting in the 
compactness range depicted in Figure [5] As before, the 
overall behaviour is in good agreement with previous 
results; the larger spreading of the various EoS in Fig- 
ure Is] is most likely due to numerical errors which are 
introduced by the coordinate system used for our simu- 
lations. 

Based on the data for nonrotating configurations, a 
linear fit leads to 



1 

fo(s) 



M3 



22.49 - 14.03 



(17) 



where again M and R are given by | [T3) . 

The number of unknown variables which determine 
frequencies and damping times of the /-mode in rela- 
tions | (10) - |[12} and ([^ - | (T7) can be reduced further by 
the well-known fact that to a very good accuracy, typi- 
cally between 5 - 7%, the Kepler-limit can be estimated 
by 



O. 



0.67 



GM 



(18) 



see e. g. f56,57|. This is only 23% larger than the Newto- 
nian value for polytropic stars. For realistic equations of 
state on the other hand, a similar empirical relation can 
be derived. In this case it was found that 



^K = C{Xs) 



GM 



with 



Xs 



2GM 
Rc^ 



and C(;t;s) =0.468 + 0.378;ts (20) 



can reproduce the original values with a relative error 
of only 1.5% ||58l|59]l. In this sense, n^ is not an in- 
dependent parameter but can be computed very accu- 
rately from the mass and the radius of the nonrotating 
neutron star 



C. Asteroseismology Examples 

The empirical relations found in this Section can be 
used in two ways. By prescribing M, R and O, one can 
easily compute frequencies und damping times of both 
the CO- and counterrotatrng mode branches for any rota- 
tion rate up the mass-shedding limit. 

On the other hand they allow to do asteroseismology, 
for example three independent measurements, two fre- 
quencies and one damping time, will lead to a robust 
estimate of mass, radius and angular frequency and will 
therefore help to restrict the range of possible equations 
of state to those in agreement with these measurements. 
Of course, this is a very idealized point of view because 
it will be very difficult to observe damping times and 
frequencies of stable oscillations; this applies both for 
the CO- and counterrotating modes. In these cases, a pos- 
sible way of estimating the feasibility of a detection is to 
set a certain threshold on the gravitational wave ampli- 
tude and relate it to the energy that has to go into the 
/-mode as it was already done in [3J for the nonrotat- 
ing case. It was shown there, that the oscillations of a 
nascent neutron star can potentially be detected from 
within the local group of galaxies or even within the 
Virgo cluster if one assumes a more optimistic estimate 
for the radiated energy. 

Before damping /growth times are included in the ex- 
amples, one might first ask what kind of information can 
be extracted by just detecting the /-mode frequencies of 
the CO- and counterrotating branch. Since the relations 
1 10 1 - ((12) and ([Ts) depend on average density and rota- 



(19) 



tion rate, one will not be able to determine mass and ra- 
dius separately but merely the combination M/R and 

n/n^. 

To give a simple example, we choose a certain back- 
ground model of the less compact EoS II sequence with 
/-mode frequencies (j\/7.n = 2.250 kHz and (Ji/1n = 
1.804 kHz in the inertial frame. Solving for the mean 
density and angular velocity yields M/R^ = 0.397 and 
n/Qj^ = 0.199 while the correct values for this par- 
ticular model are 0.467 and 0.173 respectively. As one 
can see, the relative error in average density and rota- 
tion rate is only around 15% in both cases. Of course, in 
this example mass and radius cannot be determined in- 
dependently; for this additional information about the 
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damping times need to be taken into account as it is 
shown next. 

We will address two separate questions here. The first 
one concerns the accuracy of the fittings when compared 
to the exact results and the second examines the possi- 
bility to use them for actual astroseismology, i. e. how 
accurate they constrain the neutron star parameters for 
a given set of measured frequencies and damping times. 

As first example, we choose the less compact equi- 
librium sequence of EoS P0.66 and an arbitrary value 
for the rotation rate, e. g. 0/2/1 = 0.676 kHz which 
is roughly 38% of the mass-shedding limit and corre- 
sponds to a ratio of polar to equatorial coordinate radius 
of r-n/re = 0.95. Inserting the values for mass and ra- 
dius of this particular model from Table [TV] into relation 
| (T8) leads ioVt^llnK, 1.742 kHz which is only slightly 
smaller than the correct value of Vi.^l7.n = 1.775 kHz 
computed from the background code. 

The /-mode frequencies extracted from our time evo- 
lution are given by ij\l1n = 4.192 kHz for the coro- 
tating and (JiI'^tl = 2.397 kHz for the counterrotating 
branch; both frequencies are given in the inertial frame. 
On the other hand, evaluating the fitting formulas (TO) , 
([11), |[12), ^ with the correct values for M, R and fJ 
yields ai/lTi = 4.205 kHz and 0-2/271 = 2.425 kHz re- 
spectively w^hich is an excellent match with the correct 
frequencies. 

A similar comparison can be made for the damping 
times. Using the correct /-mode eigenfunctions for com- 
puting energy and energy-loss due to gravitational radi- 
ation, i. e. relations l|5) - 0, leads to ti = 0.035 sees and 
T2 = 0.559 sees. Then again, when using the correct val- 
ues for M, R and n in relations ^ - ^ and ^ - ^ 
one arrives at ti = 0.050 sees and f2 = 0.702 sees respec- 
tively. This agrees quite well with the correct results; see 
Table IT] for a summary. 

parameter exact value value from the fit 



measurements. To continue with the previous exam- 
ple, we first look at the tuple of measurements A^i : = 
[cTi, (72, Ti) that is the two /-mode frequencies of the co- 
and counterrotating branch and the damping time of the 
corotating mode. Reversing relations ( [TO) , (TTj l and ^T6\ 
leads to a system of nonlinear equations that is solved by 
an iterative root-finding algorithm as described in [60]. 
For a reasonable choice of starting parameters, this al- 
gorithm converges with adequate precision to an esti- 
mate of M, R and Q. More specifically, from Si we get 
M = 1.44 M0, _R = 8.98 km and n/n^ = 0.386 and 
using relation (TS) with these values of mass and radius 
leads to Ci/ln = 0.669 kHz. 

One should note, that for a given tuple of measure- 
ments A4, the root-finding algorithm can in principle 
lead to other solutions as well, depending on the ini- 
tial guess for the beginning of the iterations. For exam- 
ple, if one uses (M,_R,n/nj;c) = (2.0,1.0,0.0) as start- 
ing values, the root finder converges to M = 11.54 M0, 
R = 17.98 km and n/n^ = 0.386. However, the es- 
timates for M and R are well beyond the range of ex- 
pected neutron star masses and radii so that they can 
safely be discarded although the rotation rate is matched 
perfectly. On the other hand, discarding these unphys- 
ical solutions, the nonlinear solver converges to essen- 
tially the same set of roots independently of the initial 
setup for the start of the iteration. In this example, the 
various solutions obtained with different initial guesses 
differ by less than 2%. 

Alternatively, one can also check the correspond- 
ing results when providing the measurements A^2 '■ = 
{o'i,o'i,'^i) as input data, this time with the damping 
time of the potentially CFS-unstable branch. Since now 
a different fitting function for t is used, the corre- 
sponding findings from the nonlinear root solver will 
be slightly different in general. However, for this partic- 
ular example we find practically the same values for M, 
R and O/Ox as in the first case; see Table [ill 



CLk/Iu 


1.775 
4.192 
2.397 


1.742 
4.205 
2.425 












<Ti/2n 




M 


R 


a/In 


0^/271 


a2/ln 


exact 


1.10 


8.18 


0.676 


1.775 


Tl 


0.035 


0.050 


using M-i 


1.44 


8.98 


0.669 


1.732 


T2 


0.559 


0.702 


using M2 


1.44 


8.98 


0.669 


1.732 



TABLE I. Comparison of the stellar parameters Q^, cr (both in 
kHz) and t (in sees) for a less compact EoS P0.66 model with 
O = 4.247 kHz between actual simulations and the empirical 
relations. 



TABLE II. Comparison between exact and estimated stellar pa- 
rameters from solving the inverse problem for measurements 
Ml, M2- Here, M is in units of M©, R in km and n. Ox in 
kHz. 



Overall, the accuracy by which the fitting relations 
can reproduce the exact results from actual simulations 
is reasonable. The difference is less than 2% for the var- 
ious frequencies and 25% - 40% for the damping times. 

The second issue we would like to address is the in- 
verse problem, i. e. how precise are the restrictions on 
mass, radius and angular velocity for a given triple of 



We repeated these two types of checks, i. e. compar- 
ison between exact results and fitting functions as well 
as solving the inverse problem, also for other EoS and 
larger angular velocities and a summary of these results 
for the more compact equilibrium model of EoS PI. 2 is 



given in Table III 



This time, the angular velocity of the actual model 



parameter exact value value from the fit 



nx/27r 

(72/271 
T2 



0.777 
2.259 
1.043 
0.099 
4.716 



0.806 
2.283 
1.059 
0.093 
6.650 







M 


R 


n/27r 


CiK/ln 


exact 


1.58 


15.42 


0.464 


0.777 


using M-i 


1.51 


15.33 


0.479 


0.796 


using M2 


2.43 


17.97 


0.479 


0.795 



TABLE III. Comparison between exact and estimated stellar 
parameters for an EoS PI .2 model rotating at roughly 60% of 
the mass-shedding limit. M is in units of M©, R in km and n. 
Ox in kHz while t is given in sees. 



is increased to about 60% of the mass-shedding limit. 
Consequently, the damping time of the counterrotating 
mode T2 is in the range of a few seconds already and 
it will continue to grow for more rapidly rotating mod- 
els until the CFS-instability begins to operate. Still, the 
fitting functions provided in this study can reproduce 
the exact values for the frequencies within an error of 
3% while the difference in the damping times again is 
around 6% - 40%. 

Concerning the solutions of the inverse problem, the 
two different measurements A^i and A^2 lead to differ- 
ent values for mass and radius of the neutron star in this 
example. Here, using either A^i or A42 alone over- or 
understimates somewhat the correct values for M and 
R whereas the angular velocity and the Kepler-limit is 
matched almost perfectly with any of the datasets; the 
error there is around 3%. 

Altogether we were able to demonstrate the applica- 
bility of the fitting functions provided in this study for 
a wide range of different polytropic equations of state 
and rotation rates ranging from moderate to rapid. The 
accuracy of the empirical relations decreases for rapidly 
rotating models as expected but otherwise provides re- 
sults within an error of a few percent for the frequencies. 
The dependency of the fits for the damping times on the 
mode-frequency is very strong; this is especially true for 
the potentially unstable modes and hence the estimates 
for the damping times are less accurate the closer one 
approaches the Kepler-limit. 

The solution of the inverse problem delivered rea- 
sonable values for stellar key parameters like mass, ra- 
dius and rotation rate. Typically, using the measure- 
ments Ml or A^2 alone either under- or overestimates 
the correct values for mass and radius while the rota- 
tion rate proved to be more robust to estimate by a sin- 
gle measurement alone. Ideally, the scheme presented 
here should be used iteratively in the following way: A 
first solution of the inverse problem already puts some 



constraints on the possible values for mass, radius and 
rotation rate in parameter space. Based on these con- 
straints, certain EoS and rotation rates can be discarded 
and improved fits for frequencies and damping times 
within this region in parameter space will lead to fur- 
ther restrictions for mass, radius and angular velocity. 
This procedure can then be repeated until the desired 
accuracy is reached. 

Another important issue is related to the issue of accu- 
rately determining mode parameters by an actual grav- 
itational wave detector However, a thorough discus- 
sion of measurement errors in frequencies and damping 
times is beyond the scope of this paper; we refer the in- 
terested reader to |4| and references therein. 



IV. SUMMARY 

In this work we demonstrated how one can do gravi- 
tational wave asteroseismology by using the frequencies 
and possibly the damping /growth times of the emitted 
waves from oscillating and rapidly rotating relativistic 
stars. This is possible by the empirical relations that we 
have derived and w^hich connect the frequencies and the 
damping /growth times of the oscillation modes with 
the stellar characteristics, i.e. with the mass, radius and 
rotation rate. We have actually shown that for poly- 
tropic equations of state of varying stiffness one can cre- 
ate very robust formulae connecting the observable fre- 
quency and damping times with the quantities like ro- 
tation frequency, average density and /or compactness. 

We have shown on a few examples how one can use 
the empirical formulae in order to derive the stellar pa- 
rameters. In a realistic situation when an /-mode will 
be excited, it will be possible to detect the signal at least 
from galactic sources if the mode is CFS-stable and at 
least from sources in the Virgo cluster if it is unstable 
ll6TH64l . This will be possible with the sensitivity of the 
advanced Virgo and LIGO detectors fSF, '66] and prob- 
ably even more feasible with the next generation gravi- 
tational wave telescopes such as ET (Einstein Telescope) 

IHIISZI. 

The "weak" point of the whole procedure relies in the 
approximate calculation of both frequencies and damp- 
ing times. As we already mentioned we have neglected 
the spacetime perturbations and thus there is a system- 
atic quantitative but not qualitative error in all data. 
Thus it is expected that the coefficients in relations (TO) , 
1 11 1 , I JIS] ! and 1 (16) will be affected by a proper treatment 
of the spacetime degrees of freedom. However, it is be- 
lieved that these changes will not alter the results signif- 
icantly since the relations for frequencies and damping 
times are normalized with their corresponding values in 
the non-rotating limit. These values will absorb most of 
the differences when compared with the correct results 
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in the presence of spacetime perturbations but of course, 
this is an issue that has to be addressed properly in fu- 
ture work. An additional outcome of this analysis will 
be the frequencies and damping times of the zf-modes. 

Finally, the empirical relations found in this study 
have been derived for polytropes of varying stiffness 
which are able to mimic the global properties of realis- 
tic equations of state; see for example EoS A and EoS II 
which are polytropic fits to tabulated EoS. Realistic hot 
equations of state are the best candidates for newly born 
neutron stars and have a higher chance of becoming un- 
stable prior mutual friction completely suppresses any 
instability |68| but currently our code is unable to per- 
form time-evolutions of rapidly rotating neutron stars 
for generic tabulated data. Preliminary studies regard- 
ing this issue are promising but in a very early stage. 
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Appendix A: Numerical Procedure 

This work is a continuation of our previous efforts 
||43l |44| and relies on the foundations laid therein. We 
will therefore briefly summarize the crucial parts of the 
previous studies that are needed for the computation of 
the damping times here. 

We are numerically solving the relativistic hydrody- 
namics equations, linearized around background equi- 
librium configurations of uniformly rotating neutron 
stars. For this purpose, the time-evolution of the fluid 
perturbations is performed in a cylindrical coordinate 
frame {q,l,,(p) which is comoving with the neutron star, 
surface-fitted for all rotation rates and where the metric 
takes the form 



ds2 = e-2^ 



e^^ (dq^ + df\ + W^dcp- 
e™{dt + adcpf. 



(Al) 



Here, the metric potentials U, k, W and a depend on q 
and I, only and are obtained by solving the generalized 
TOV-equations for axisymmetric equilibrium configu- 
rations and a perfect-fluid energy-momentum-tensor 
For computational purposes, the physical domain of the 
simulations 



V = [{q,l,),cp = const] 



(A2) 



is mapped onto a rectangular grid 

r=[s,t)e [0,1] X [0,2], 

{q = q{s,t),l = as,t))eV], (A3) 

where the barotropic fluid equations for the perturbed 
velocity and pressure are discretized and, together with 
appropriate boundary conditions, are numerically in- 
tergrated by using an Iterated Crank-Nicholson scheme 
with an additional amount of Kreiss-Oliger dissipation. 
Actually, we do not directly evolve the fluid perturba- 
tions but certain combinations of hydrodynamical and 
metric variables. This considerably reduces the com- 
plexity of the differential equations and simplifies the 
boundary treatment of the computational domain. More 
specifically, for azimuthal mode number m our time- 
evolution variables are given by 



f^{q,^,t)e'''^V = {e + p)]Ne^5u^ 
f2{q,^,t)e'"-'f = ie + p)We^3u^ 
Mq,^,t)e""'P = ie + p)Su^ 
H{q,C,t)e'""l' =c^e^Se, 



(A4) 



where /i, /2, /a, H are integrated in time and Su^, Su^^, 
Suq,, Se are the perturbed fluid velocities and energy- 
density respectively. Furthermore, p is the background 
pressure, e the unperturbed energy-density and Cg the 
speed of sound w^hich can be computed analytically 
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from the equilibrium configuration for polytropic equa- 
tions of state. 

After a successful time-integration, which is typically 
cancelled after 50 - 70 ms, one can use Fast Fourier 
Transforms to extract oscillation frequencies and eigen- 
functions of any mode one wants to study. This data 
is then used for further post-processing; for example in 
order to compute damping times. 

For this, one has to adapt equations dSl - |I7|| to the 
computational domain T. The corresponding relations 



are 



E = 



1 



/J 



pSu'Sul + ^Sp* 



\detj{s,t)\dsdtd(l) 

(A5) 



for the energy contained within a mode. 



dE 
It 



= -ai{cri + mn)N24\^D22(' 



(A6) 



for the energy-loss due to gravitational radiation and 

D22 = / y/7Spr^Y^2\detJ{s,t)\dsdtdcl) (A7) 

for the mass-quadrupole moment. Here 

7= {W^exp{-2U)-a^exp{2U))exp{-AU + Ak) 

(A8) 
is the determinant of the spatial 3-metric from l(Al|, 
J{s,t) is the Jacobian matrix of the coordinate mapping 
| |A3) and r is the radial distance. Due to the azimuthal 
decomposition of the perturbation variables, the inte- 
gration in ^-direction is trivial and raising the covariant 
fluid perturbations with the inverse metric yields 



SU^SUn 



Su^Sur 



{e + . 



)Wexp{k) 
1 



(e + p)Wexp(A:) 



1/1 



1/2 



(A9) 



Su^Su 



f 



expjU) 



I/3I 



Note, that energy-density and pressure vanish at the 
stellar surface, the same applies for the metric potential 
W along the rotation axis. These critical points might 
lead to numerical problems when evaluating the expres- 
sions in | |A9) . However, in practice it turns out that 
our time-integration ensures the correct behaviour of 
the perturbation variables in order for equations I A9 1 to 
remain finite. 

Furthermore, the adiabatic condition leads to 



Sp 



H 



exp(U) 



(AlO) 



for the pressure perturbation, whereas the correspond- 
ing change in density can be computed once an equation 



of state is specified. In this study, we are considering 
polytropic equations of state (EoS) which take the form 

p = Xp^+i/~ where £ = p + Np. (All) 

Here K is the polytropic constant, N the polytropic ex- 
ponent and F = 1 + 1 / N the polytropic index. In grav- 
itational units {G = c = Mc-j = 1), K^^^ can be used as 
scaling factor and in this new system, one arrives at 



Sp = -p^/^-^Sp. 



(A12) 



Finally, the Jacobian matrix needs to be computed. This 
can either be done by first interpolating Q{s,t), l,{s,t) 
onto the computational domain T, followed by a finite- 
difference routine or alternatively by interpolating the 
pseudo-spectral expressions for the transformation co- 
efficients fl^y := {dxi/dyi^; Xj E (s, t),i/j- E ((), ^)}and 
using the inverse function theorem to obtain the Jaco- 
bian of the inverse function. Due to interpolation errors 
in the first approach, the subsequent finite-difference 
scheme leads to highly non-smooth results for the Jaco- 
bian / and the error introduced there can be as high as 
30%; we are therefore favouring the second approach. 

Eventually, standard finite integration schemes are 
deployed to evaluate the integrals in the expressions for 
energy and energy-loss | |A6| and I A7l. 

As an example and in order to get an impression of 
how the extracted eigenfunctions are arranged on the 
computational domain. Figure [6] shows the power spec- 
tral density of the scalar perturbation variable H, see 
| |A4) for the definition. 

There, the value of s labels the radial coordinate which 
starts at the origin of the star for s = and terminates 
at the stellar surface for s = 1 . The use of surface-fitted 
coordinates ensure that the surface is always located at 
s = 1, even for rapidly rotating models. On the other 
hand, t acts as an angular coordinate ranging from f = 
at the rotation axis above the equatorial plane (i. e. ^ > 
0) to f = 1 at the equatorial plane and to f = 2 at the 
rotation axis below the equatorial plane (^ < 0). 

It is then clear, that the eigenfunction depicted in Fig- 
ure |6] can indeed be identified with the fundamental 
I = m = 2 oscillation mode. It has no nodes in the ra- 
dial direction and the angular pattern matches the scalar 
spherical harmonic Y22 ^ sin 9 which has a maximum 
for 9 = n/2, i.e. the equatorial plane, as well as global 
minima along the rotation axis. The eigenfunction also 
vanishes at the origin and along the surface of the star; 
this is due to boundary conditions and the special choice 
of the time-evolution variables, see [44J for more details. 

One should also keep in mind, that the power spec- 
tral density essentially measures the energy contained 
within a mode; this is at least true for the velocity per- 
turbations /i, /2, /s. Any auxiliary normalization coeffi- 
cients that are included from the Fast Fourier Transform- 
algorithms cancel out in the final equation ||4| for the 
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FIG. 6. Power spectral density of H (normalized units) for the 
quadrupolar /-mode and a nonrotating background configu- 
ration. 



damping time. This means, that we can directly use the 
data output of the eigenfunction extraction routine for 
computing the energy and energy-loss of the fundamen- 
tal mode. 



Appendix B: Background Configurations 

In this study we treat neutron stars as perfect-fluid 
objects that obey a polytropic equation of state | |A11) . 
Naturally this is a rather crude approximation, neglect- 
ing a variety of micro-physical effects such as the true 
internal constitution of neutron stars (i. e. the distribu- 
tion of baryons, leptons, optionally also hyperons and 
kaons) or the influence of a finite temperature as well as 
superfluidity and the existence of a solid crust (which 
will become important for temperatures around 10^*^ K), 
see [69 To]. We also do not account for the influence of 
magnetic fields which affect oscillation modes only for 
very high field strengths, see e. g. llTll . 

In this sense, a simple polytropic equation of state 
parametrizes our ignorance about the true microphys- 
ical description of neutron star matter w^hich is still un- 
known presently. However, one can nevertheless use 
relativistic polytropes to cover the wide range of ex- 
pected neutron star masses and radii Il72l [73l and this is 
our proposed strategy here. We utilize a variety of poly- 
tropic EoS which are used in large parts in several other 
simulations of neutron star oscillations, either perturba- 
tively or non-linear II14I 139^461 174 1, two of them (EoS A 
and EoS II) actually are polytropic fits to tabulated data. 
Excluding this last two equations of state, the naming 
convention EoS P# has been chosen here, w^here # is re- 
placed by the polytropic index N. 



For each equation of state, we choose two different 
models; one rather close to the maximum allowed mass 
and another less compact one. A brief summary of ba- 
sic stellar parameters for the non-rotating configurations 
can be found in Table ITVl 



EoS 


K 


r 


M (Ms) 


Re (km) 


11 


1186 


2.34 


1.91 
1.50 


11.62 
13.19 


A 


1528 


2.46 


1.61 
1.25 


9.47 
10.54 


P0.66 


1000 


5/2 


1.29 
1.10 


7.60 
8.18 


Pl.O 


100 


2.0 


1.51 
1.40 


13.32 
14.07 


P1.2 


35 


11/6 


1.58 
1.32 


15.42 
18.38 


P1.4 


20 


12/7 


1.91 
1.83 


20.69 
23.17 



TABLE IV. EoS parameters and basic stellar properties for 
the nonrotating models. K is given in dimensionless units 
(G = c = Mq = 1), M is the gravitational mass and R^ the 
equatorial circumferential radius. 

Starting from the nonrotating configuration and by 
keeping the central rest-mass density fixed, the ratio 
of polar to equatorial coordinate radius rp/ff is suc- 
cessively decreased, typically by a factor of 0.05, until 
the mass-shedding limit is reached. In addition, sup- 
plementary equilibrium models were constructed espe- 
cially for resolving the regime of slow rotation since for 
example moving from rp/fg = 1.0 to fp/rg = 0.95 for 
EoS P0.66 already means to reach roughly 40% of the 
Kepler-limit. 

Since the centrifugal force now supports the pressure 
in sustaining gravity, rotating models attain much larger 
masses. In Table |V] we summarize fundamental stel- 
lar parameters for the maximally rotating equilibrium 
models. 

As it can be inferred from the data in this Table, the 
less compact models of any sequence can be deformed 
to a greater extent than the more compact ones as one 
would expect. This effect is of course more pronounced 
for softer equations of state. 

In order of increasing stiffness, the equations of state 
investigated in this paper can be ordered P0.66 < A < 
II < Pl.O < P1.2 < P1.4. This can also be seen from 
the mass-radius relations depicted in Figure IT] where in 
addition the uniformly rotating background configura- 
tions are referenced as well. 

Typically, the increase in mass is not as pronounced 
as the increase in the radius when turning on rotation. 
While the percental increase in the gravitational mass 
can be as low as around 10% in our sequences of con- 
stant rest-mass density, the corresponding increase in 
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EoS 


rp/r. 


a^lin 


M (Mo) 


Rt, (km) 


n 


0.564 


1.393 


2.61 (+37%) 


15.64 (+35%) 




0.546 


0.998 


1.94 (+29%) 


18.77 (+42%) 


A 


0.558 


1.759 


1.97 (+22%) 


12.71 (+34%) 




0.537 


1.282 


1.64 (+31%) 


15.05 (+43%) 


P0.66 


0.555 


2.180 


1.58 (+22%) 


10.24 (+35%) 




0.540 


1.775 


1.42 (+29%) 


11.50 (+41%) 


Pl.O 


0.579 


0.966 


1.80 (+19%) 


18.57 (+39%) 




0.570 


0.853 


1.70 (+21%) 


19.83 (+41%) 


P1.2 


0.597 


0.777 


1.81 (+15%) 


21.39 (+39%) 




0.593 


0.543 


1.57 (+19%) 


26.19 (+42%) 


P1.4 


0.62 


0.538 


2.13 (+12%) 


28.55 (+38%) 




0.61 


0.445 


2.07 (+13%) 


32.64 (+41%) 



TABLE V. Basic stellar properties of the maximally rotating 
configuration and percental increase compared to the nonro- 
tating case. rplVf, is the ratio of polar to equatorial coordinate 
radius and CiY^lln represents the Kepler-limit, given in kHz. 
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FIG. 7. Mass-Radius diagram for the six EoS and in our sam- 
ple. The twelve different background sequences (two for every 
EoS) are marked as circles which branch off their correspond- 
ing non-rotating equilibrium curve. 



equatorial circumferential radius is always larger than 
at least 30%. 



Appendix C: Consistency Check 

Here, we compare our results for the damping time 
of nonrotating polytropes with literature values in full 
general relativity. 

In [j75l/ the authors compute damping times and fun- 
damental mode frequencies for several relativistic poly- 
tropes with different polytropic indices. In Newtonian 
theory, damping times and mode frequencies for poly- 
tropes only depend on mass and radius of the neutron 
star and on the polytropic index N of the particular 



equation of state but not on the polytropic constant K 
or the central energy-density This scale invariance is no 
longer valid in general relativity but it can be shown, 
that the relativistic calculations match with the Newto- 
nian ones in the limit of small masses and large radii. 

More specifically, if one defines a dimensionless pa- 
rameter c-c via 



R 



fGM\ 



(CI) 



then the relativistic results for a F = 2-polytrope ap- 
proach the scale-invariant Newtonian value of c^ = 8.46 
in the low-mass-limit. 

For this purpose, a sequence of non-rotating back- 
ground models was computed in 1'751 with N = 1 and 
X = 100 km^ and central rest-mass densities ranging 
from Q, = 3x 10^^ - 0.05 x lO^^ g/cm^. Then, damp- 
ing times and oscillation frequencies of the fundamental 
quadrupolar mode were computed and compared with 
their Newtonian counterparts. 

We repeated the calculation with our code and an ad- 
justed equilibrium sequence of EoS Pl.O which has the 
same polytropic index as the background models but a 
different value of K. Adopting K = 100 km'^ leads to 
more compact neutron stars with a maximum mass of 
^max = 1-1 Mq and Rmax = 7.2 km; compare with Fig- 
ure U\ One should also keep in mind that the results 
in [75 1 are obtained by solving a complex eigenvalue 
problem. Here, we perform the calculation with our 
method of computing the eigenfunction from a time- 
evolution of the fluid perturbations and evaluating the 
integrals | |A5) , | |A6) . These are two completely different 
approaches and an agreement between these two meth- 
ods would be a strong indication for the accuracy of our 
procedure. Figure l8] shows a comparison between the 
values of Cj obtained by these two approaches. 
The overall agreement between the different methods is 
very good. Using a second-order fit for our results leads 
to a value of c-^ w 8.95 in the Newtonian limit compared 
to Cr = 8.46 found in [75J. Especially for models with 
low M/R, the least compact object in our simulations 
has a mass of M = 0.034 M© and a radius oi R = 12.5 
km, the /-mode frequencies are in the range of several 
hundred Hertz which impairs a proper eigenfunction 
extraction with our code. 

The last remark concerns the limit of large masses in 
Figure l8] Although we also use a polytropic equation 
of state, both are implemented slightly different here 
and in |75|. Here, a relation of the form | |A11 1 is con- 
sidered, w^hich properly describes an ideal gas under- 
going adiabatic processes, while the authors in 1 75 1 use 
p = Kq^^^'^ with e = (? in the relativistic case. This de- 
scription of the fluid neglects the pressure contribution 
to the energy density and permits the speed of sound 
to become potentially larger than c for all values of the 
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FIG. 8. Comparison between literature values and our results 
for c-z (see JCTJ for the definition). In the limit of small com- 
pactness, the Newtonian value of c-^ = 8.46 is approached. 
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polytropic index N, see also [76]. Both prescriptions lead 
to the same stellar models in the limit of small compact- 
ness, because the energy density is dominated by the 
rest mass in this case and the contribution of the pres- 
sure is negligible. In more relativistic cases, the pres- 
sure contribution is noticeable and leads to a decrease 
in the mass of the most massive star which is d5n-iami- 
cally stable to radial oscillations. This is the reason why 
our sequence of equilibrium configurations already ter- 
minates at Mmax ^ 1-1 Mq while the sequence in Il75l 
can reach up to Mmax ~ 1-3 Mq. However, this has only 
a modest effect on the computation of c^- 
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